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ABSTRACT 

This paper contains details of computer models of shower 
development which have been used to investigate the 
experimental data on shower cores observed in the Leeds 
35m 2 and Sacramento Peak (New Mexico) 20m 2 arrays of 
current-limited spark ("discharge' 1 ) chambers. The 
simulations include predictions for primaries ranging 
from protons to iron nuclei (with heavy nuclei treated 
using both superposition and fragmentation models) . 


1. Introduction 

Analysis of photographs of 10 l4 -10 16 eV air shower cores in arrays 
of current-limited spark chambers at Leeds (1020 gem" 2 ) and Sacramento 
Peak ( 730 g enf 2 ) has provided data on particle distributions within a 
few metres of the shower axis in the form of density spectra and density 
averages measured at known fixed distances from the shower centre (Ash, 
1985a, Hodson et al., 1983a, b, 1985). Comparison of these results with 
Monte Carlo simulation predictions by Ash (1985a, b) has demonstrated the 
failure of the scaling model. Details of the calculations are presented 
in this paper (see also Ash, 1984) . 

2. Hadron collision models 

Simulation techniques based in the laboratory frame were used, and 
the models generated only pions and leading nucleons; nuclear target 
effects were neglected. The energy-splitting algorithms of Hillas (1981) 
were used to generate "radial scaling" secondary energy spectra in pion 
collisions (Fig. la) and nucleon collisions (Fig. lb),- except that for the 
latter, the flat leading nucleon spectrum given was replaced by a more 
accurate form (Hillas, 1979a): the elasticity X - Enuc/ E o was sampled 

using X = 0*76a + 0*243** (instead of X = a) where a is a random number 
(0*0. . . . 1*0) . (If, as a result, E nuc was less than the nucleon rest 
energy me 2 then E nuc = me 2 ) . For E 0 » me 2 the resulting average 
elasticity is equal to 0*42. 

Transverse momenta were sampled from the energy-independent 
distribution h(pt) dpt “ pte~ pt ' p ° with p Q equal to 0*24 GeV/c for 
secondary pions (in both pion and nucleon collisions) and equal to 
0*28 GeV/c for leading nucleons. (If pt ^ p then the secondary was 
considered to be going backwards and was dropped from the cascade.)- 

Mean free paths for nucleon- 1 air 1 nucleus collisions (Fig. 2) based 
on nucleon-nucleon cross sections rising as &n 2 s at high energy were as 
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(a) Pion collisions (b) Nucleon collisions 
Fig. 1 Particle production spectra in simulated 100 GeV hadron-'air' 
nucleus collisions. X is the secondary energy divided by the 
projectile energy, both measured in the laboratory frame 


given by Hillas (1979b); those for pion collisions were obtained by 
scaling up the nucleon mean free paths by (114/86) at all energies. 


3. Nuclear fragmentation models 


The superposition model" was u 
number A, energy Eprim. was taken as 
showers with energy E prim /A. A more 



Fig. 2 Inelastic pp cross- 
sections and corresponding 
mean-f ree-paths in air 
used in simulations 


sed i.e. a shower with primary mass 
a co-axial superposition of A proton 
realistic treatment of primaries 
ranging from helium nuclei to iron 
nuclei was given in the one- 
dimensional M f ragmentation model" 
'which is based mainly on emulsion 
data and follows quite closely the 
Bartol model (Gaisser et al., 1982) 
with which it is in good agreement 
for the average distribution in the 
atmosphere of the "first"' inelastic 
collisions of primary nucleons. 
Excellent agreement with Freier and 
Waddington T s (1975) simulation 
predictions of the average rate of 
release of nucleons from projectile 
nuclei as a function of atmospheric 
depth is obtained using the 
Cleghorn nucleus-nucleus cross 
sections, adopted by them, instead 
of the larger cross-sections of 
Gaisser et al. (1982) normally used. 


4. Propagation model 

A non-isothermal atmosphere was used. Coulomb scattering and 
ionisation loss (2.2 MeV per g cm -2 ) were applied to hadrons, but were 
neglected for projectile nuclei in the "fragmentation model". 
Geomagnetic effects were neglected. 
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5, Electromagnetic cascade function 

An average development function based on unpublished Monte Carlo 
calculations for photon-initiated cascades by A.M. Hillas was used. The 
electron density at 1020 g cm“^ (Leeds) for cascades due to vertical 
photons was 


o(r) 


1 dn(r) 
2*r dr 

and 


(1 ♦ r/45)~ 3 ra“ 2 with 

a - 0* 460/ ( 1 ♦ 0 • 380/ x) 
f - l - 0- 72(1 - y/t) 


, . a t(l ♦ tnf/O-72) 

n(r) - — c 

/y 

X - E.r/5-8093 

y - tn ( l ♦ 2 * 390x) 


E is the photon energy in GeV, and t is the vertical thickness (r.l.) 
between the photon production point (it 0 decay) and the observation level. 
A radiation length of 36*1 gcnT 2 was used. If a cascade was inclined this 
was taken into account by multiplying the density calculated from the 
formula by the cosine of the zenith angle and taking t equal to the slant 
thickness, and r equal to the perpendicular distance from the cascade 
axis . 

The cascade function has been used directly for each photon (EMC 
Model I) or in conjunction with a simple Monte Carlo technique to simu- 
late fluctuations in the formation of the first e + e“ pair (EMC Model II). 

Mountain level calculations involved rescaling the above electro- 
magnetic cascade lateral distributions according to the lower air 
density. 


6. Simulation of the electron distribution in the shower core 

(a) Calculations for Leeds 

A map of the core region was obtained by superposing electromagnetic 
cascades, taking into account the direction and location of their axes. 

In each shower simulation the shower primary was directed at the centre 
of a square of side 10* 2m located in the observation plane. The square 
was subdivided into (51) 2 equal squares, and the counts of electrons 
deposited in each of them determined using the cascade function. 

Computing time consumed in doing this was saved by ranking photons so as 
to concentrate the most detailed calculations on photon cascades with 
the densest "cores", and by "thin sampling" (Hillas, 1981) of low energy 
hadrons. Calculations for photon energies down to 1 GeV from tt° decay 
could therefore be made on a routine basis. 

After calculating the electron counts the simulation program 
searched for a shower centre (centre of symmetry) in the central (7*6) m 
of the mapping area (most of the shower centres, located to within ^ 5 cm, 
fall within 0.5m of the centre of the mapping area). The shower was 
discarded if a shower centre was not found in the "search area", and also 
if the electron density, averaged over 2x2 adjacent square elements of 
the mapping area, did not exceed 50 m -2 somewhere in the mapping area. 

If the shower was not discarded then densities at various radial 
distances from the shower centre were obtained, as in the experimental 
work (Hodson et al.» 1983a, b, 1985) , from counts made in an annular 
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grid", in this case, "superposed" over the mapping area. The grid used 
consisted, as in the experiment, of a series of concentric circles, 
starting with a radius of 0 25 m and then with radii increasing in 0*5 m 
steps. For the above experimental work the annuli were divided into 
equal areas of rc/16m 2 . However for the simulations it was more 
convenient to divide the annuli angularly each into 16 sections 
("sampling bins") by common radii from the centre. The difference 
between the two procedures is not significant: both give densities for 

the same radial widths. 

Counts were obtained in the "sampling bins" with the "annular grid" 
centred on the shower centre. The central count, calculated by integra- 
ting the electron density over the area of the central bin, gave the 
density corresponding to the experimental central density p(0). Counts 
in the annuli at 1*0, 2*5 and 4*0 m were derived from the mapping area 
counts and gave densities corresponding to the experimental densities 
P(1*0), P(2- 5) and p(4»0) respectively. These results were supplemented 
by counts, determined as above, in the remaining annuli out to 4*5m and 
by densities calculated at points in the mapping plane further out. The 
radial density distribution was thus sampled at intervals of 0*5 m out 
to 7-0 m from the shower centre. 

(b) Calculations for Sacramento Peak 

Both the calculations and the experimental analysis are made using 
the corresponding Leeds- level grids with dimensions in array space scaled 
up according to the lower air density found at mountain level. 

7. Construction of density spectra etc. 

Predictions of the average lateral distribution and density spectra 
at fixed distances from the shower centre have been obtained for a wide 
choice of primary composition and energy spectra. A technique of 
weighting the simulated showers according to the primary fluxes assumed 
was used to obtain results for different primary energy spectra. 
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